Wormhole formation in dissolving fractures 



P. Szymczak 

Institute of Theoretical Physics, Warsaw University, Hoia 69, 00-618, Warsaw, Poland 

A. J. C. Ladd 

Chemical Engineering Department, University of Florida, Gainesville, FL 32611-6005, USJ^\ 

We investigate the dissolution of artificial fractures with three-dimensional, pore-scale numeri- 
cal simulations. The fluid velocity in the fracture space was determined from a lattice-Boltzmann 
method, and a stochastic solver was used for the transport of dissolved species. Numerical simu- 
lations were used to study conditions under which long conduits (wormholes) form in an initially 
■ rough but spatially homogeneous fracture. The effects of flow rate, mineral dissolution rate and 

^ ' geometrical properties of the fracture were investigated, and the optimal conditions for wormhole 

formation determined. 

0\ '. I. INTRODUCTION 

i-C ■ A number of experimental and numerical studies of dissolution in fractured or porous rock have established that the 
evolving topography of the pore space depends strongly on the fluid flow and mineral dissolution rates. Remarkably, 
there exists a parameter range in which positive feedback between fluid transport and mineral dissolution leads to the 
spontaneous formation of pronounced channels, frequently referred to as "wormholes". Spontaneous channeling of a 
^p. reactive front has been shown to be important for a number of geophysical processes, such as diagenesis @, [H[ j melt 
c/2 ' migration 0, l20l .l56l . Wf\ , terra rosa formation [63| , development of limestone caves [H, Hl[ , and sinkhole formation by 
salt dissolution [761 ] . Further details can be found in review articles and books on reactive transport and geochemical 
self-organization, [e.g. [II H BE III- 

Wormholes play an important role in a number of geochemical applications, most notably CO2 sequestration IIP . 
10], [55[ , risk assessment of contaminant migration in groundwater [3J] and stimulation of petroleum reservoirs [30l l49j |. 
Selecting the optimal flow rate is an important issue in reservoir stimulation, so as to achieve the maximum increase in 
permeability for a given amount of reactant [Til. l33l l35l. IHol. |72| . If the acid is injected too slowly, significant dissolution 
occurs only at the inlet, and the permeability of the system remains almost unchanged. At the other extreme of high 
injection velocities dissolution tends to be uniform throughout the sample. However, the increase in permeability 
is again insignificant, since the reactant is consumed more or less uniformly throughout the fracture, making only 
an incremental change to the permeability. Moreover, some of the reactant may escape unused. The most efficient 
y—{ • stimulation is obtained for intermediate injection rates, where the reactive flow self-organizes into a small number of 
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distinct channels, while the rest of the medium is effectively bypassed. This focusing mechanism leads to much more 



efficient use of reactant, since the development of channels causes a large increase in permeability with a relatively 



' small consumption of reactant. 



Experimental studies of wormhole formation have used a variety of porous systems; plaster dissolved by water 



18 1, limestone cores treated with hydrochloric acid 44;] and salt-packs dissolved with under-saturated salt solution 
56 j . Recently, a variety of dissolution patterns in single rock fractures have been reported [2l], [2^, [25|, [2^, [36|, [7. 



depending on the chemical and physical characteristics of the fracture-fluid system. The physico-chemical mechanisms 



behind the pattern formation are not yet understood in detail. Linear stability analysis has been used to investigate 
the conditions required for the break up of a planar dissolution front [111 . 143 . l69( , but these results only pertain to 
the initial stages of channel formation, where the front perturbations are small. The later stages of channel evolution 
are strongly nonlinear and here numerical methods are needed. The numerical models used to study wormholing 
in porous media fall into four broad categories. (1) Single wormhole models 0, |48| . in which a channel of a given 
shape is created in the porous matrix. The reactant concentration field inside the wormhole is determined and the 
growth velocity calculated. (2) Darcy-scale models [H, [zH based on continuum equations with effective variables 
such as dispersion coefficients, Darcy velocity and bulk reactant concentrations. These average variables replace the 
microscopic diffusion constant, fluid velocity, and reactant concentration. (3) Network models [H,[44|], which model 
fluid flow and dissolution in a network of interconnected pipes, where the diameter of each network segment or pipe 
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is increased in proportion to the local reactant consumption. (4) Pore-scale numerical simulations 0, HH, HH, [HI • In 
these calculations the equations for fluid flow, reactant transport and chemical kinetics are solved in an explicitly 
three-dimensional pore space. Although computationally intensive such models provide detailed information on the 
evolution of the fluid velocity, reactant concentration and topography without invoking effective parameters such 
as mass transfer coefficients. This is the approach followed in the present work, however in the context of fracture 
dissolution. 

In studies of fracture dissolution, and particularly in theoretical investigations of cave formation, a one-dimensional 
model of a single fracture is frequently used [e.g. Q, [2(| [27], H3|- The fracture aperture (the distance between the 
rock surfaces) is assumed to depend on a single spatial variable, the distance from the inlet. Although analytically 
tractable, one-dimensional models cannot account for wormhole formation, and thus they are only relevant at the 
extremes of high and low flow rate where the dissolution is expected to be uniform. In two-dimensional models of 
dissolving fractures [HI, SI| , the fluid velocity and reactant concentration are averaged over the aperture of the 
fracture. The key simplifications are the Reynolds (or lubrication) approximation for the fluid velocity [l[ and the use 
of effective reaction rates. These models are technically similar to Darcy-scale models, with the local permeability 
determined by the aperture; they produce realistic-looking erosion patterns and correlate positively with experimental 
results [22j. However the Reynolds approximation may significantly overestimate the flow rate 8, 64, 68], especially for 
fractures of high roughness and small apertures. Moreover, under certain geological and hydrological conditions, large 
pore-scale concentration gradients develop and in such cases volume averaging can introduce significant errors [59l . l60j , 
sometimes not even capturing the correct reaction direction. 

The most fundamental approach is to directly solve equations for fluid flow, reactant transport, and chemical 
kinetics within the fracture space. This approach was pioneered by Bekri et al. [f|, who solved the flow and transport 
equations using finite-difference schemes. Better resolution is offere d by lattice Boltzmann methods, which have been 
used in dissolution simulations at the pore scale 0, H HI IH HI S M, HI and at the Darcy-scale [H, H3]- Here 
we combine velocity field calculations from an implicit lattice-Boltzmann method [85| with a transport solver based 
on random walk algorithms that incorporates the chemical kinetics at the solid surfaces [HI]. Advances in numerical 
algorithms for flow and transport allow us to simulate systems of relevance to laboratory experiments without resorting 
to semi-empirical approximations [831 ]. In the simulations reported here, as many as 50 interacting wormholcs have 
been studied (c.f. Fig. [17)) , comparable to systems modeled by state of the art Darcy-scale simulations [3], while 
maintaining pore-scale resolution. 

We have investigated wormhole formation in a simple artificial geometry, where one of the fracture surfaces is 
initially flat, and the other is textured with several thousand randomly placed obstacles. The geometry is similar to 
that studied experimentally by Detwiler et al. [12, HH and has the advantage that it shows no discernible long-range 
spatial order. The correlation length is of the order of the distance between the obstacles, and is much smaller than the 
system size. Although it lacks the self similarity of natural fractures, it provides a useful starting point for numerical 
analysis of wormholing, since the initial structure contains no nascent channels. 

The aim of this paper is to discover the range of conditions under which long conduits form in a initially rough 
but spatially homogeneous fracture. Dissolution was studied under conditions corresponding to a constant pressure 
drop across the sample and to a constant flow rate. Constant pressure drop is representative of the early stages of 
karstification [26|, whereas constant flow rate is more relevant for reservoir stimulation (30j . In this context, we have 
numerically determined the conditions needed to maximize the permeability increase for a given amount of reactant. 

We investigate the effects of flow rate (characterized by the Peclet number) , mineral dissolution rate (characterized 
by the Damkohler number), and geometrical properties of the fracture. The Peclet number measures the relative 
magnitude of convective and diffusive transport of the solute, 

Pe = vh/D, (I) 

where v is a mean fluid velocity, h is the mean aperture and D is the solute diffusion coefficient. In this work, 
v = Q/(Wh) is related to the volumetric flow rate, Q, and the the mean cross sectional area, Wh, where W is the 
width of the fracture. The Damkohler number, 

Da = k/v, (2) 

relates the surface reaction rate to the mean fluid velocity. The relevant geometric characteristics are harder to 
quantify. Hanna and Rajaram [4l| argued that the key geometrical factor determining the intensity of wormholing is 
the statistical variance of the aperture field, cr, relative to the mean aperture, / = cr/h. Here, we present numerical 
evidence that the total extent of contacts between the surfaces may play an important role as well. These contacts 
need not be load bearing; the dynamics remains qualitatively the same if the surfaces are sufficiently close that the 
local fluid flow is strongly hindered. The transition between uniform dissolution and channeling seems to occur rather 
sharply in our simulations, at the point where the contact regions make up 5% — 10% of the total fracture area. 



3 



inlet 

manifold 




FIG. 1: The geometry of the experiment: a corrugated glass surface (upper) is matched with a soluble flat plate (lower). The 
plates are held in a fixed position and reactant flows from an inlet manifold designed to produce a uniform concentration and 
flow field at the inlet. 



II. NUMERICAL MODEL 



To investigate channel growth and interaction in a dissolving fracture, we use a pore-scale numerical model [83j in 
which the fracture space is defined by two-dimensional height profiles h u (x, y) and h (x, y) representing the upper and 
lower fracture surfaces. The velocity field in the fracture space is calculated by an implicit lattice-Boltzmann technique 
[85| , while the transport of dissolved species is modeled by a random walk algorithm, which efficiently incorporates 
the chemical kinetics at the solid surfaces 82] . The fracture surfaces are discretized into pixels and the height of each 
pixel is eroded in response to contacts by tracer particles; for the results reported below either 200 x 400, 400 x 400, 
or 800 x 800 pixels were used. The time evolution of the velocity field and the local aperture variation in the fracture 
are determined by removing small amounts of material at each step, and recalculating the flow field and reactant 
fluxes for the updated topography. 



A. Flow field calculation 

In laboratory-scale fractures, the Reynolds number is less than 1 [23], [25, 28]; it is also small during the initial stages 
of cave formation [13, (zH- Thus, inertia can reasonably be neglected, and fluid motion is then governed by the Stokes 
equations 

V ■ v = 0; ??V 2 v = Vp, (3) 

where v is the fluid velocity, r\ is the viscosity and p is the pressure. The velocity field in the fracture has been calculated 
using the lattice-Boltzmann method with "continuous bounce-back" rules applied at the solid- fluid boundaries [S6| . 
The accuracy of these boundary conditions is insensitive to the position of the interface with respect to the lattice, 
which allows the solid surface to be resolved on length scales less than a grid spacing; thus the fracture surfaces 
erode smoothly. It has been shown [87| that the flow fields in rough fractures can be calculated with one-half to one- 
quarter the linear resolution of the "bounce-back" boundary condition, leading to an order of magnitude reduction in 
memory and computation time. A further order of magnitude saving in computation time can be achieved by a direct 
solution of the time- independent lattice-Boltzmann model [S5| . rather than by time stepping. These improvements 
have previously allowed us to calculate velocity fields in laboratory-scale fractures [HI] . The calculation of the flow 
field in a 200 x 400 fracture takes about 1 minute at the beginning of the dissolution process, and up to 15 minutes 
during the final stages of dissolution. The processor was a single core of an Intel Pentium P4D clocked at 3GHz. 



B. Solute transport modeling 

Solute transport in the fracture is modeled by a random walk algorithm that takes explicit account of chemical 
reactions at the pore surfaces. The concentration field is represented by a distribution of tracer particles, each 
representing n solute molecules. We use a standard stochastic solution to the convection-diffusion equation, [SB], [13] 

d t c + v • Vc = £>V 2 c. (4) 

in which individual particles are tracked in space and time, 

n(t + St) = Ti(t) + v(ri)St + V2DStT. (5) 

The flow field, v(r), is derived from the implicit lattice-Boltzmann simulation and T is a Gaussian random variable 
of zero mean and unit variance. The fluid velocity at the particle position is interpolated from the surrounding grid 
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points, while the time step St is chosen such that the displacement in one step is smaller than O.lSx, where Sx is 
the grid spacing. To account for chemical erosion at the fracture surfaces, we calculate the dissolution flux at each 
boundary pixel, assuming a first-order surface reaction 

Jl = k(c s - c ), (6) 

where c s is the saturation concentration, cq is the local concentration at the surface and k is the surface reaction rate. 
The solute flux is normal to the surface, and forms a boundary condition to the transport solver, 

-D(V±c)\ = nk(c s -co), (7) 

where n points into the fluid, and Vj_ = nn • V. The notation V(. . .)|o indicates a gradient at the surface. 

It is convenient to introduce a reactant concentration field C, so as to simplify the boundary condition in Eq. ([7]); 
in the present context 

C = c s -c (8) 

is the undersaturation, measuring the deviation of c from the saturation concentration. The boundary condition ([7J 
is then 



D(V±C)\ = nkC (9) 

In a different situation, for instance dissolution of a fracture by a strong acid, the reactant field, C, is the acid 
concentration. It is most convenient to define C so that the boundary condition always takes the form of Eq. ; the 
convection-diffusion equation ((4]) is the same in both cases. 

The drawback of the classical random walk method Q is that a very large number of particles must be tracked 
simultaneously, so that the concentration near the rock surface can be determined accurately enough to obtain a 
statistically meaningful dissolution flux. However, there is a considerable simplification in the case of linear dissolution 
kinetics, where it is possible to derive a single-particle stochastic propagator that satisfies the boundary condition in 
Eq. §§§ [82| . The diffusive part of the particle displacement in the direction perpendicular to the fracture surface (z) 
is then sampled from the distribution 

G d (z,z',5t) = 1 ( p -(z-zT/lDSt +p -(*+zT/lDSt\ 

V^DSt V / 
_ fe e k(kSt +z+z >)/n Er{c ( z + * + M 8t\ 

In Eq. (jTUJ) , z' and z are the distances of the tracer particle from the surface at the beginning and at the end of the 
time step respectively. The boundary condition ^ implies that the amount of reactant represented by a single tracer 
particle n(t) decreases in time according to 



n{t + 5t)=n{t) G d {z,z',St)dz. (11) 
Jo 

The integral can be calculated analytically, 

n(t + 5t)/nit) = e^'+^Erfc (^=^f) (12) 



+ Erf 



z 



VlDSiJ ' 

and for k > the amount of material represented by the tracer is reduced, 

n(t + 5t)/n(t) < 1. (13) 

In order to apply Eq. (|10[) to a complex topography, the timestep 5t must be limited, so that within each step a particle 
only samples a small portion of the fracture surface, which then appears planar. In our simulations V DSt ~ 10~ 2 /ioj 
where is the initial mean fracture aperture. 

The side walls of the fracture (parallel to the flow direction) are solid and inert; thus reflecting boundary conditions 
for the solute transport (Eq. [7] with fc = 0) are imposed at y = and y = W. At the fracture inlet (x = 0), a 
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reservoir boundary condition of constant concentration C — C{ n is applied, while a saturation condition C = C ou t = 
is assumed at the outlet boundary, x = L. These boundary conditions are implemented according to the algorithms 
described in [8l|, which contain a number of subtleties. Because mineral concentrations in the solid phase are typically 
much larger than reactant concentrations in the aqueous phase, there is a large time-scale separation between the 
relaxation of the concentration field and the evolution of the fracture topography. We therefore make a quasi-static 
approximation, solving for the time-independent velocity and concentration fields in a fixed fracture geometry. The 
quasi-static approximation may break down in cases where the reactant is much more concentrated and the reaction 
kinetics are fast; acid erosion by HC1 is a possible example of this. 

The steady dissolution flux in the fracture can be calculated by tracking individual tracers, using the following 
algorithm [si!, [12]: 

1. Sample the initial position of a tracer particle within the inlet manifold indicated in Fig. [TJ Assign the initial 
number of reactive molecules represented by a tracer, 

n(0) = -^C in : (14) 

Vo is the volume of the inlet manifold, Ntot is the total number of particles to be sampled, and Cm is the inlet 
concentration of reactant. 

2. Propagate the particle for a single time step St, according to Eq. ([5|). If it comes within a cutoff value z c , of 
any surface element, then sample the diffusive part of the particle displacement perpendicular to the wall from 
Eq. (fT0|) . change n(t) according to Eq. ([13]), and increment the dissolution flux counter at the surface element 
(i) closest to the particle by 

AJj = n(t + ft)-n(«) 
bidt 

where Si is the area of the surface element. 

3. The random walk described by repeating step 2 many times is terminated when n(t)/n(0) falls below a preset 
threshold, or when the particle leaves the system through the inlet or outlet. Random walks are also terminated 
when the particle fails to enter the fracture at the first step, but these must be counted, even though they do 
not contribute to the erosion flux. 

4. Upon completion of N to t random walks, remove material from the fracture walls in proportion to the accumulated 
fluxes, Ji. The total amount of material is chosen to be sufficiently small that the evolution of fracture topography 
appears continuous. 

5. The cutoff distance was set to z c — \Q\j2DSt, N to t was of the order of 10 6 , and the threshold below which a 
tracer is deleted was 10 _6 n(0). 

The above scheme can be used to calculate concentration profiles in large fractures, since it is more computationally 
efficient than a typical stochastic algorithm where the local concentration field is needed to determine the dissolution 
flux. In that case the number of tracer particles (Ntot) required for statistically significant erosion rates is several 
orders of magnitude larger. 

The time evolution of the velocity field and local aperture are determined by iteration, removing small amounts of 
material at each step. The dissolution-induced aperture change in the fracture over the time, Ah, is related to the 
mean dissolution flux J — J^. JiSi/J2i &i by 

A -^JA^ (ig) 

Csol " a q 

where c so i is the concentration of the solid component and v aq , v so i are the stoichiometric numbers of the aqueous 
and solid species. It is more computationally efficient to keep Ah constant in each erosion cycle, and then increment 
the time accordingly, 

At= 2f^^L. (17) 

J V so l 

In the simulations reported here, Ah = 0.0lh n . Test calculations with smaller values of Ah (down to Ah = 0.001/io) 
confirmed that the patterns are insensitive to the magnitude of the erosion step in this parameter range. 
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We will use a dimensionless timescale, based on the time for high flow rate dissolution in a parallel channel. In 
this idealized system, the reactant concentration is everywhere uniform and equal to the inlet concentration Cj„. The 
spacing between the plates is set to ho, the initial value of the mean fracture aperture, and the reference reaction 
rate is chosen so that the product of Peclet and Damkohlcr numbers is unity, PeDa = kho/D = 1. We define the 
characteristic time r for a plate to erode by ho- 

_ _ % CsqI v aq 
D Ci n v sol 

Equation (|17[) can then be rewritten in terms of a dimensionless time At/r, 

At DC,„ Ah 



t h 2 J 



(19) 



It is important to stress that the above model contains no free parameters or effective mass-transfer coefficients. 
Instead the fundamental equations for fluid flow, reactant transport, and chemical kinetics are solved directly. The 
simulations incorporate the explicit topography of the pore space, and the transport coefficients-viscosity, diffusivity, 
and reaction rate-are determined independently. 



C. Validation of the numerical model 



The numerical model has been validated [83] by comparison with experimental data obtained with an identical 
initial topography (23|. The experimental system was created by mating a 99 x 152mm plate of textured glass (spatial 
correlation length of ~ 0.8mm) with a flat, transparent plate of potassium-dihydrogen-phosphate (KDP). The relative 
position of the two surfaces was fixed during the experiment, eliminating the effects of confining pressure, which are 
hard to control experimentally 28] and even harder to model numerically [87| . The fracture was dissolved by an 
inflowing solution of KDP at 5% undersaturation. High spatial resolution data (1192 x 1837, 0.083 x 0.083mm pixels) 
was obtained for the evolution of the local fracture aperture as a function of spatial position [231 ] . The experiments 
were conducted at two different hydraulic gradients, corresponding to initial mean velocities v = 0.029 cm s _1 and 
0.116 cm s~ x . The other parameters characterizing the system are the diffusion coefficient of KDP in water, D = 
6.8 x 10 -6 cm 2 s _1 , the initial mean aperture, ho = 0.0126 cm, and the reaction rate, k — 5.2 x 10~ 4 cm s" 1 . Note that 
in [HI the reaction rate was erroneously reported to be smaller by a factor of two, k — 2.6 x 10~ 4 cm s _1 ; however 
the correct value was used in all of the calculations reported there. The Peclet and Damkohler numbers calculated 
for these parameters are, from Eqs. Q and (J2J), Pe = 54, Da = 0.018 (v = 0.029 cm s" 1 ) and Pe = 216, Da = 0.0045 
(v = 0.116 cm s -1 ). 

Figure [5] shows a comparison of the dissolution patterns obtained by simulation and experiment; the initial topogra- 
phies in the simulation and the experiment were the same. The contour levels in the figure represent the change in 
height of the lower (dissolving) surface (Fig. [1} as a function of time 

Ah(x, y;t) = h l (x,y,0)- h l (x, y, t). (20) 

At higher flow rates, unsaturated fluid penetrates deep inside the fracture and dissolution tends to be uniform 
throughout the sample (right panels in Fig. [2]), while at the lower flow rate erosion is slower and inhomogeneous 
(middle and left panels in Fig. \2§ ■ The dissolution front is unstable to fingering [691 ] , since an increase in permeability 
within a channel enhances solute transport through it, reinforcing its growth. As dissolution proceeds, the channels 
compete for the flow and the growth of the shorter channels eventually ceases. At the end of the experiment, the flow 
is focused in a few main channels, while most of the pore space is bypassed. 

The experimental and numerical dissolution patterns are similar. At low Peclet number, the dominant channels 
(Fig. [2]) develop at the same locations in the simulation and experiment, despite the strongly nonlinear nature of the 
dissolution front instability. While there are differences in the length of the channels, relatively small changes (of 
the order of 10%) in the diffusion constant, D, or rate constant, k, can lead to comparable differences in the erosion 
patterns. Our results suggest that the simulations are capturing the effects of the complex topography of the pore 
space; a more extensive and quantitative discussion, including histograms of aperture distributions at different Peclet 
numbers, can be found in [83j . 



FIG. 2: (Color online) Erosion of the lower surface (initially flat) during dissolution of a laboratory-scale fracture. Dissolution 
patterns for Pe — 54, Da = 0.018 are shown at Ah = ho/2 in the left panels and at Ah — ho in the middle panels; the right 
panels show dissolution patterns at Pe — 216, Da = 0.0045, Ah = ho. The simulations are shown in the lower panels and the 
corresponding experimental results are shown in the upper panels. The lightest/red shading indicates the deepest erosion while 
dark gray/blue indicates the least erosion; the intermediate colors are yellow (higher) and green (lower). The flow direction is 
from left to right. 



The computational model described in Sec. Ill Bl was used to simulate dissolution in artificial fractures, with numeri- 
cally generated topographies. Initially, the lower surface of the fracture is flat, while the upper surface is textured with 
several thousand identical cubical protrusions (parallelepipeds of height 2Sx and base 35x x 3Sx, where Sx is the pixel 
size). The protrusions were placed on a square lattice and then randomly shifted by ±Sx in the lateral (y) direction, 
which eliminates all the straight flow paths between the inlet and outlet. The resulting fracture has short range 
spatial correlations, and no discernable long-range structure. The fracture geometry can be characterized statistically 
by the fractional coverage of protrusions, £. If the obstacles span the entire height of the fracture aperture, the initial 
geometry has a relative roughness 



where a is the variance of the aperture field. Most of the simulations discussed here have been carried out for £ = 0.5, 
which corresponds to a relatively rough fracture / = 1; however we also investigated smoother fractures, with Q as 
small as 0.025 (see Sec.|X|. 



A typical initial geometry (£ = 0.5) is shown in the left panel of Fig. [3J the integrated (two-dimensional) velocity 
field, 



is shown in the center panel. Reactive fluid enters from left side and exits from the right, while no-slip boundaries 
are imposed on the other surfaces. The setup and initial topography resemble the experiment in Detwiler et al. [23[, 
but here we allow both surfaces to dissolve, which speeds up dissolution and channel formation. In this case the 
time-dependent erosion depth is defined in terms of the combined change in height of both upper and lower surfaces, 



III. ARTIFICIAL FRACTURE GEOMETRIES 




(21) 




(22) 



Ah(x,y;t) = [ 



h u (x,y,t)-h l {x,y,t)] 
h u (x,y,0)-h l (x,y,0)] . 



(23) 



Initially there are no discernible channels (center panel) and the velocity field shows only short-range spatial corre- 
lations. During dissolution, large scale variations develop from small fluctuations in the initial porosity. In the final 
stages of dissolution (right panel in Fig. [3]) the channeling is very distinct, but the size of these spontaneously formed 
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FIG. 3: (Color online) Initial distribution of obstacles (dark pixels) in the artificial fracture (left); initial flow field (center); and 
the flow field after an increase in mean aperture equal to its initial value, Ah = ho (right). The flow field, v 2d = \/v 2d ■ v 2d , 
Eq. I)22[l. is averaged over the local aperture. The lightest/red shading indicates the highest velocity, while dark gray /blue 
indicates the lowest; the intermediate colors are yellow (higher) and green (lower). The flow direction is from left to right. 



channels are not related to the initial pore size distribution, which is highly uniform. The growth in mean aperture, 

A ^ (t) = iw j J A M*,»>*)dyd*, ( 24 ) 

will sometimes be used as a measure of elapsed time; it is then normalized by the initial mean aperture h . 



IV. CHANNELING AS A FUNCTION OF PECLET AND DAMKOHLER NUMBER 



Figure 2] illustrates typical dissolution patterns at a mean erosion depth Ah = 2ho, over a range of different Peclet 
(Pe) and Damkohler (Da) numbers. Changes in erosion patterns map more uniformly to variations in the inverse 
Damkohler number, Da^ 1 , than variations in Da. We therefore use Da -1 as the independent variable in most of our 
plots, although, in conformity with normal practice, we discuss the results in terms of variations in Da. 

For small Pe and large Da the reactant saturates (C = 0) near the injection face. After a fast initial dissolution 
of material at the fluid inlet, the reaction front propagates extremely slowly, as there is almost no unsaturated fluid 
penetrating inside the fracture. On the other hand, when the reaction rate is sufficiently slow (Da < 1/100), or the 
flow rate sufficiently high (Pe > 500), unsaturated fluid penetrates deep inside the fracture and the whole sample 
dissolves almost uniformly. Channeling is observed for moderate values of Peclet and Damkohler numbers, Pe ~ 10 
and Da > 1/100. Here nonlinear feedback plays a decisive role. A perturbation in the reaction rate at the dissolution 
front increases (for example) the local permeability, which in turn increases solute transport and therefore the local 
dissolution rate. The increasing flow rate reinforces the initial perturbation and the front becomes unstable, developing 
pronounced channels where the majority of the flow is focused, while most of the pore space is eventually bypassed. 
The interaction between channels, important in the later stages of dissolution, is discussed in Sec. ILXl 

Above a threshold Damkohler number, Da > 1, the erosion patterns are largely determined by the Peclet number. 
Additional data (not shown) demonstrates that there is no difference in dissolution patterns at Da = 1, Da = 10 
and Da — > oo, except at small Peclet numbers, Pe < 1. A similar range of Damkohler number (0.1 < Da < oo), 
has been reported in other numerical studies [79( as a regime where the length of a single dissolving channel in a 
two-dimensional porous medium becomes independent of reaction rate; we will return to this point in Sec |Vj In the 
mass-transfer-limited case, where the surface reaction rate is high enough that the overall dissolution process no longer 
depends on reaction rate, the stochastic modeling may be simplified by imposing an absorbing boundary condition, 
C = 0, at the fracture surfaces, corresponding to the limiting case Da — > oo; these results are shown in the leftmost 
columns of Fig. [U 

For smaller Damkohler numbers, Da < 1, the dissolution patterns become dependent on both Pe and Da. The 
Peclet number controls the number of channels, with the spacing between them decreasing with increasing Pe. On the 
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FIG. 4: (Color online) Erosion of the lower surface (initially flat) at constant pressure drop (left) and constant injection 
rate (right), for different Peclet and Damkohler numbers. The lightest/red shading indicates the deepest erosion while dark 
gray/blue indicates the least erosion; the intermediate colors are yellow (higher) and green (lower). The flow direction is from 
left to right. 



other hand, for fixed Pe, a decrease in Damkohler number results in a more diffuse boundary between the channels 
and the surrounding porous matrix. When the Peclet number is less than one, diffusive transport becomes more 
important than convection, even in the flow direction. In this regime, the dissolution patterns are determined by the 
product of Peclet and Damkohler number, 



(25) 
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FIG. 5: (Color online) Dissolution patterns in the diffusive regime, for Pe =1/2 (upper row) and Pe = (lower row). The 
values of (PeDa)" 1 = D/kho are marked. 
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FIG. 6: Phase diagram describing characteristic dissolution patterns as a function of Peclet and Damkohler number. The 
points mark the values of Peclet and Damkohler number corresponding to the patterns shown in Fig. [4] and the crosses mark 
the KDP fracture experiments shown in Fig. [2] Since the scales are logarithmic, the points corresponding to Da^ 1 — are 
marked at Da -1 — 0.1, to which they are similar. The arrows indicate the direction the points move in the phase diagram 
during dissolution at constant pressure drop (solid arrow) and constant flow rate (dashed arrow). 



which gives the relative magnitude of reactive and diffusive fluxes. Figure [5] compares the dissolution patterns 
for Pe = 1/2 with those for the purely diffusive case, Pe = 0. The differences are rather slight, except at very 
small reaction rates; the reaction front propagates slowly and stably inside the fracture, with the penetration length 
decreasing with increasing PeDa. 

These results are summarized in a phase diagram, Fig. O where the values of Peclet and Damkohler number 
corresponding to the patterns shown in Fig. g] are marked, together with the points corresponding to the KDP 
fracture experiments by [23| (see Sec. Ill C|) . Although the geometry of the KDP fracture is different from the obstacle 
fracture geometries considered here, the dissolution patterns captured at comparable Pe and Da are nevertheless 
similar. For example, the KDP dissolution pattern at Pe — 54, Da — 0.018 (Fig. [2|) may be compared with the 
artificial fracture at Pe = 32, Da = 0.025 (Fig. gj). 

A feature of reactive flows, as compared with other pattern forming systems, is that the dimensionless numbers 
characterizing the flow and reaction rates are changing throughout the course of the dissolution [19| . Indeed, in 
constant pressure drop simulations (left panel of Fig. 2]) , both the total flow and the mean aperture change during 
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the course of dissolution; thus the point in the phase diagram representing the initial system moves towards larger 
Pe and smaller Da values as the dissolution progresses (solid arrow in Fig.[6|). 

Constant pressure drop conditions are representative of many groundwater flow systems, including the early stages 
of karstification [26j |. However, in a number of technological applications, e.g. acidization of petroleum reservoirs [30l |. 
the control variable is the injection rate of reactive fluid, Q. In that case, the Peclet number remains constant 
throughout the dissolution process, since 

vTi Q 

Pe =D=WD- (26) 
On the other hand, the Damkohler number then increases in proportion to h (dashed arrow in Fig. [5]), since 

k khW 

Da = -_= ——. 27 

v Q 

The results of the constant injection rate simulations are presented in the right panel of Fig. 0] While the differences 
are not dramatic, channels formed in constant pressure drop simulations are noticeably more diffuse at their tips. 
This is consistent with the analysis in Sec. IVIi where it is shown that the thickness of the dissolution front near 
the channel tip is proportional to Da^ 1 . The Damkohler number decreases in the course of constant pressure drop 
simulations while it increases during constant flow rate runs, which leads to different front thicknesses at the tip. 
This observation agrees with laboratory experiments [44| , where it was observed that wormholes in acidized limestone 
formed by constant pressure drop dissolution become more highly branched at later times than those formed at 
constant flow rate. 

We have verified that the results shown in Fig. 2] are independent of the length of the fracture domain. In Fig. 
the dissolution pattern in a 200 x 400 pixel fracture are compared with a longer 400 x 400 domain at the same 
Peclet and Damkohler numbers, Pe = 8, Da = 1/10. The initial topographies of both systems were identical in 
the 200 x 400 pixel inlet region, and the comparison was made when the same volume of material had been eroded 
from each sample. We find that the dissolution patterns in both systems are very similar, up to the point where the 
dominant channels reach the outflow of the smaller domain, as can be seen in Fig. [Jj In Sec. HXi we describe how 
the simple hierarchical growth pattern of the competing channel system then becomes disrupted, with large pressure 
gradients developing at the tips of the leading channels, causing them to split (see also Fig. [TTjl . Here we simply wish 
to point out that, in the wormholing regime, the results shown in Fig. [5] are insensitive to further increases in the 
length of the domain. However, in the uniform-dissolution regime, the overall length of the fracture is important. The 
average undersaturation decays exponentially along the flow direction and so will eventually reach a region where the 
solution is saturated and no dissolution occurs. In this sense it seems that the distinction between uniform dissolution 
and surface inundation is merely a matter of scale. When the length of the fracture is comparable to the depth of 
penetration, dissolution appears uniform, but if the fracture is much longer, then dissolution is limited to a small 
region (relative to the length) near the inlet. 
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V. EFFECTIVE REACTION RATE 



The phase diagram in Fig. [6] can be better interpreted in terms of an effective reaction rate, which takes account 
of the interplay between mass transfer and chemical kinetics. In a two-dimensional description of dissolution [22l |. 
the fracture is locally approximated by two parallel plates separated by a distance h(x,y,t). In the reaction-limited 
regime, diffusive timescales are much shorter than reactive ones, kh/D <C 1, and the concentration field is almost 
uniform across the aperture. The reactant flux is then given by 

J = kC w kC (28) 

where C(x, y) is the cup-averaged concentration. In the opposite limit, kh/D 3> 1, dissolution becomes mass-transfer 
limited and an absorbing boundary condition, Cq = 0, may be assumed at the fracture walls. The dissolution flux, 



is determined by the mass transfer coefficient kt, 



J = ktC, (29) 



k t = Sh%-, (30) 



where dh is the hydraulic diameter of the system and Sh is the Sherwood number. For parallel plates dh = 2h and 
Sh = 7.54 [|. 

In the general case when kh/D ~ 1, the reactant flux may expressed in terms of C by equating the dissolution flux 
(1281) to the mass-transfer flux ll29l) 



J = kC = k t {C - Co) (31) 

Solving for Co in terms of C gives 

J = k ef fC (32) 

with 

keff = (33) 

There are two approximations made here. First, the Sherwood number (|30[) depends, in general, on the reaction 
rate, k. However, the variation in Sh is relatively small [13)113, bounded by two asymptotic limits: constant flux, 
where Sh — 8.24 for parallel plates, and constant concentration, where Sh = 7.54. Here we take the approximate value 
Sh = 8 in order to estimate Da e ff. Second, we have neglected entrance effects, which otherwise make the Sherwood 
number dependent on the distance from the inlet, x. However, the entrance length (defined as the distance at which 
the Sherwood number attains a value within 5% of the asymptotic one) is negligibly small, L x w 0.008dhPe [2!|, at 
least in the early stages of the dissolution. Expressions for the effective reaction rate coefficient analogous to (13"3"|) have 
been proposed previously 2^, 27 , [72, [74[ , but a somewhat different approach was employed by Howard and Groves 
|4(| and Hanna and Rajaram 41 1 where the smaller of the two rates k and k t was used for k e f f . 

The effective reaction rate, k e ff (Eq.[331 can be used to construct an effective Damkohler number, Da e ff — k e ff/v, 

2Pe 

Da- e j f = Da- 1 + —, (34) 

which remains finite even when the reaction rate becomes very large, Da — > oo; contours of constant Da e ff are shown 
in Fig. [8] The phase diagram of dissolution patterns in the Pe — Da e ff plane, shown in Fig. [9j is simpler than the 
corresponding phase diagram in Pe — Da (Fig.[S]). In particular, uniform dissolution can now be uniquely associated 
with small Da e ff, whereas in Pe — Da variables it corresponds to either small Da or large Pe. The introduction 
of Da e ff also explains the independence of the dissolution patterns on the microscopic Damkohler number in the 
Pe > 1, Da > 1 regime discussed in Sec. lIVI At higher Peclet numbers a large change in Da corresponds to a relatively 
small change in Da e ff. For example, at Pe = 32, Da e ff changes from 0.125 when Da — * oo to 0.111 when Da = 1. 

The general features of the phase diagram agree with experimental and numerical studies of wormhole formation in 
quasi-two-dimensional porous media |35J . However, in these simulations the transitions between different dissolution 
patterns corresponded to fixed values of either Peclet or Damkohler number; in other words, the boundaries in the 
Pe — Da e ff phase diagram were perpendicular to the axes. In our case the line between surface inundation and the 
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FIG. 9: Phase diagram describing characteristic dissolution patterns as a function of Peclet and effective Damkohler number. 

wormhole regime is not horizontal; at low Peclet numbers, the critical value of Peclet number at which channels are 
spontaneously formed decreases with decreasing Da e f f . The discrepancies at low Pe may be caused by a transition to 
a more three-dimensional flow in our fracture simulations. Althou gh t he initial geometry shares many features with 
the two-dimensional porous media considered in other studies [35], 72], the third dimension plays a more significant 
role as the dissolution proceeds, particularly in a large Da e ff, small Pe regime, where the penetration of the reactive 
fluid is very small. In this regime, the solutional widening of the fracture at the inlet can be more than one order of 
magnitude larger than the mean aperture growth in the system and thus the system ceases to be quasi two dimensional. 

VI. FRONT THICKNESS 

A detailed examination of Fig. 2] reveals a number of qualitative features of the developing channels. At high 
Damkohler numbers and low Peclet numbers the channels are very distinct, with sharp well-formed boundaries 
between dissolved and undissolved material. As the Peclet number increases the channels become more diffuse, but 
only in the flow direction. The lateral thickness of the channels is almost unchanged, while along the flow direction 
a sharp transition is replaced by a gradual change in dissolution depth, which takes place over almost the whole 
channel length; this is especially pronounced above Pe — 30 in the constant pressure-drop case and Pe — 100 in the 
constant flow rate case. The second qualitative feature is that at smaller Damkohler numbers, the channels become 
laterally diffuse as well, as can be seen from the broad blue regions at the dissolution front when Da < 0.1. Insight 
into the characteristics of channel formation can be gained from a simple model based on a Darcy-scale description 
of a wormhole. 

Consider the tip of a channel containing reactive fluid, with depth-averaged concentration C, entering the undis- 
solved medium. At the leading edge of the wormhole, the Darcy-scale equation for reactant transport takes the 
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FIG. 10: Sketch indicating the characteristic dimensions of a wormhole front. The parameters \ x and \ y indicate the extent 
of the dissolution into the porous matrix. 



form [6l|,l79|, 

9C_ d*C dC 2k eff 

at -° x dx* V dx h ' ( j 

where D x = D(l+ f3 x Pe) is the dispersion coefficient along the flow direction, which includes the effects of fluctuations 
in the fluid velocity through the coefficient fi x ~ 0.5 (72J. The last term describes the concentration loss due to disso- 
lution in the medium ahead of the reaction front, assuming the fracture here is undissolved and may be approximated 
by two parallel surfaces separated by ho (Sec. |V|). Then the erosion flux at both upper and lower surfaces is k e ffC 
and the rate of change in concentration is —2k e f fC /ho- The stationary solution of ([55]) is [79| . 

C{x) = C tip e- X > x * (36) 

where Cu p is the reactant concentration at the tip of the channel. 

The characteristic thickness of the dissolution front, A^, is given by (rill. |79| 



, _ 2D X 



8D x k eff ^ 1/2 



h v 

Along the flow direction, convective effects are typically much stronger than diffusive ones, and Eq. (|3"T|) simplifies 

hp ho ( 1 2Pe\ 



(37) 



2Da e ff 2 \Da k>n J 

where X x is the thickness of the front ahead of the wormhole tip (see Fig. [TO]). Convective effects are small in directions 
transverse to the flow, and setting v = in the analogue of Eq. (|35[) . gives the transverse thickness of the reaction 
front, 

,, = (MoV /2 = So (l±MiV /J , (39) 



2 fee// / \2PeDa e f f 

where the dispersion coefficient f3 y ~ 0.1 j72j. 

In the reaction-limited regime (PeDa <C 1), ~ ho/2Da, independent of Peclet number. This scaling prediction 
can be observed in Fig. [4] most clearly at Da — 1/40, for both constant pressure drop and constant flow rate 
conditions. The extent of the front roughly corresponds to the blue (color figure) or dark grey regions in front of the 
wormhole. At Da =1/40 these diffuse regions extend roughly 20-30% of the length of the channel (~ 40 — 60ho), 
independent of Pe. The scaling with Da is approximately linear, which can be seen most clearly when comparing 
Da = 1/40 with Da = 1/10 at constant pressure drop; at constant flow rate the extent of the dissolution front is too 
small to measure at Da = 1/10. In the mass-transfer-limited regime (PeDa ^> 1), « hoPe/Sh. The extent of 
the reaction front at Da = oo (Fig. [5]) depends roughly linearly on Pe, in agreement with the simple scaling law. A 
quantitative comparison is suspect for several reasons: the difficulty in defining and measuring the extent of the front, 
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FIG. 11: (Color online) Examples of channels of a similar length at different Peclet numbers (from left to right) Pe — 2, 
Pe = 8, and Pe — 32. At the two higher Peclet numbers we used Da — oo, but at Pe = 2 we used a lower Damkohler number, 
Da = 0.1, since the channels do not form at Da = oo in this case. 



changes in Peclet and Damkohler numbers as dissolution proceeds, and local variations in fluid velocity. Nevertheless, 
the Darcy-scale model qualitatively explains the key results of the simulations shown in Fig [4] 

The behavior of the transverse thickness of the front is more complicated, because of the dispersion term in Eq. (|39[) . 

— 1/2 

In the reaction-limited regime, X y = ho [(1 + (3 y Pe) / (2PeDa)] 1 , and the line PeDa = 1 is the division between 
sharply defined and laterally extended wormholes (see Fig. 0J. At sufficiently high Peclet numbers the dispersion 
term should eliminate the dependence of X y on Pe, but the value of PeDa is then too large for an observable lateral 
extension of the reaction front. An exception may be the case Pe = 128, Da = 1/160 with constant flow-rate 
conditions (right panel of Fig [4]). This is a transition case between wormholing and uniform erosion, and here we can 
see a significant lateral spreading of the front. Much more clearly defined is the growth in lateral thickness as PeDa 
gets smaller. This again can be most readily seen at Da = 1/40; here the lateral extension of the front gets more 

— 1 /2 — 

pronounced as Pe is reduced. Finally, in the mass-transfer-limited case, X y = h [(1 + l3 y Pe)/ShY^ < h , and the 
reaction front is always limited to the region ahead of the tip, as is observed at Da = oo (Fig. |4]). 
Panga et al. [72j argue that the ratio of these two length scales 

A = ^ (40) 
X x 

determines the aspect ratio of the wormhole, with the strongest wormholing predicted to occur when Awl. The 
latter criteria matches quite well with the patterns observed in the simulations. However, examination of Fig.|4]shows 
that the aspect ratio of the wormholes does not correlate well with A. For example, at constant Pe, the wormhole 
diameter is only weakly dependent on Da (Fig. 0]), and the aspect ratio remains more or less constant, whereas the 
front thickness increases considerably with decreasing Da. Thus the aspect ratio of the reaction front, as measured 
by A, does not necessarily control the aspect ratio of the wormhole itself. The dependence of the wormhole shape 
on Peclet number is analyzed from a different perspective in the next section, based on a microscopic mass balance 
within the wormhole. 



VII. CHANNEL SHAPE 



The aspect ratio of the dissolving channels increases with increasing Peclet number as can be seen in Fig. 1111 For 
a fixed length, the typical wormhole diameter decreases approximately as Per 1 / 2 . This scaling can be understood 
by noting [79j], that, in the mass-transfer-limited regime, the shape of the channel depends on the interplay between 
diffusive transport normal to the flow, and convective transport along the flow direction. The channel boundary is 
parameterized by the curve R(x), where R is the distance of a point on the boundary from the channel center line; the 
geometry is illustrated in Fig.[TJ] The depleted concentration at the point {x, R{x)} diffuses towards the center of the 
channel with a timescale that can be estimated by solving a two-dimensional diffusion equation for the concentration 
C(r, t), inside a circle with an outer boundary condition C(R,t) = 0. This gives an asymptotic (in time) dependence 
of the concentration at the center line 



c(o,t) ~ e - a « Dt / R \ 



(41) 



16 



R(x) 




< ► 

L 

FIG. 12: The radius of a wormhole, R, as a function of the longitudinal coordinate, x. The arrows mark the characteristic time 
scales of diffusive and convective transport. 



where ao ~ 2.4 is the first zero of the Bessel function Jo- Thus the characteristic time for depleted reactant to diffuse 
from the channel boundary R(x), is given by 

The diffusing concentration field is simultaneously advected by vAt to the tip of the wormhole located at {L, 0} 
(Fig. [12]). Equating the two timescales, 



L — x R(x 



,2 



At = ^27T> ( 43 ) 

leads to an expression for the channel shape 



a%D(L-x)^ 1/2 



R{x) = p-^-^ j = Roy /T^JZ ) (44 ) 
where 

R»=(^)" 2 (45) 

is the radius of the wormhole (Fig. [T2")) . 

The analysis in this section uses a microscopic model of reactant transport within the wormhole, in contrast to the 
previous section (Sec. IVII[) where we considered reactant transport in the porous matrix. Thus, in the flow direction 
we assume that convection is dominant, whereas in the transverse directions the transport is molecular diffusion rather 
than dispersion. This is based on a picture of a wormhole as a region of low (or vanishing) porosity, so that the flow 
is roughly parabolic, constrained by the upper and lower fracture surfaces and by the boundary of the wormhole; the 
fluid velocity in the wormhole is much higher than in the surrounding matrix. Along the flow direction, diffusion is 
enhanced by Taylor dispersion, 

£>H =D(l+p T Pe 2 ), (46) 

where j3r is a coefficient that depends on geometry, but is ~ 0.005. At moderate Peclet numbers, Pe > 10, Taylor 
dispersion is comparable to or larger than molecular diffusion, but even at the highest Peclet numbers in these studies 
(Pe ~ 100) dispersion makes a negligible contribution to the axial transport of reactant. The timescale for convective 
transport over the wormhole length L is still much smaller than the diffusive timescale, even if Taylor dispersion is 
included, 
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FIG. 13: The dimensionless parameter, F = Rq (Pe/Zio^) , characterizing the wormhole shape, as a function of a Peclet 
number, Eq. 



The last inequality is valid up to at least Pe ~ 10 2 for channel lengths in excess of 10/io- 

In this simple wormhole model, Eq. (|4l))) . the radius of the channel scales like Pe -1 / 2 , while its shape is parabolic. 
To connect the theory with the numerical simulations, we have calculated the dimensionless quantity, 

T = R (Pe/h L) 1/ \ (48) 

which Eq. (|45[) predicts will have a universal value of do rs 2.4. In making these comparisons we must take into account 
that the analysis leading to Eq. (|45|) considers only individual channels. However, channel competition is an essential 
component of the overall dynamics [84j , in which the longer channels drain flow from the shorter ones, limiting their 
growth. A detailed analysis of channel competition is given in Sec. IIXI here we aim to limit the effects of channel 
competition by focusing on the longest channels, which remain active throughout the dissolution process. Fig. 1131 
shows the dimensionless ratio T, Eq. (|48|) . for dissolving fractures in the mass transfer limit. In accordance with the 
above remarks, only the three longest channels in each fracture were measured. The numerical results confirm that 
r is nearly universal, independent of Peclet number and the choice of channel. Moreover, the numerical values of T 
are close to 2.4, but this may be accidental given the limited precision of the numerical data and the simplicity of the 
model. 

In the opposite case, when the process is reaction- limited and Pe 3> 1, the reactant concentration in the entire 
wormhole is nearly uniform and equal to the inlet concentration, Ci„. Interestingly, in this limit, the shape of the 
wormhole is also parabolic [65| . 



VIII. PERMEABILITY EVOLUTION AND OPTIMAL INJECTION RATES 



The evolving permeability in a dissolving fracture can be defined by the relation 

Q = _KWh* 1 (49) 
M 

where K is the permeability, /i is the viscosity and Who is the initial cross-sectional area of the fracture. Figure [14] 
shows the permeability as a function of time for dissolution at constant pressure drop. At low Peclet numbers (Pe = 2 
in Fig. I14p . the flow rate through the sample only increases significantly at very small Damkohler numbers, where the 
dissolution is uniform throughout the fracture. At higher Damkohler numbers, surface inundation occurs and the flow 
rate remains nearly constant throughout the simulation. The behavior changes dramatically as the Peclet number is 
increased. At Pe > 20, the dependence of the flow rate on Damkohler number is reversed: now the flow increases 
most rapidly for larger values of Da. This is because wormhole formation, triggered in this parameter range, becomes 
amplified as the reaction rate increases. 

A characteristic feature of systems that exhibit channeling is the rapid growth of the flow rate when the dominant 
channels break through to the outflow end of the fracture; breakthrough is indicated by the near vertical lines in 
Fig. [T3J At moderate Peclet numbers, Pe 10, uniform dissolution and wormholing compete with each other, as 
can be seen from the permeability evolution at Peclet number Pe = 8 (Fig. [T4|). The permeability at first increases 
fastest at the lowest Damkohler number Da = 1/160, since the unsaturated fluid penetrates deeper inside the sample. 
However, once channels begin to form, the permeability in the higher Damkohler systems increases very rapidly, and 
eventually overtakes the Da — 1/160 system. 
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FIG. 14: Permeability as a function of time during dissolution at constant pressure drop; results are shown for Pe — 2, 8, 32, 
and 128. The lines correspond to: Da = 1/160 (dotted), Da = 1/40 (dash-and-dotted), Da = 1/10 (dashed), and Da = oo 
(solid). 

At constant flow rate, shown in Fig. Qj3 the permeability increases more slowly overall than at constant pressure 
drop, where dissolutional opening of the fracture is enhanced by the increasing flow rate (note the different time scales 
in Figs. [T4l and H5|) . Positive feedback is particularly strong near breakthrough, which manifests itself in steeper K{t) 
curves in the constant pressure drop simulations. Another difference between the constant pressure drop and constant 
flow rate can be observed at Pe — 8. At constant pressure drop the permeability growth is faster for Da — ► oo than 
for Da = 0.1, whereas for constant flow rate this order is reversed. This is again connected with the increasing flow 
rate in the course of dissolution at constant pressure drop. Since Pe = 8, Da — » oo lies near the borderline between 
wormholing and surface inundation regimes, even small differences in flow have a pronounced impact on the speed at 
which channels propagate. 

A quantitative description of channeling under constant flow conditions is also important in carbonate reservoir 
stimulation, where the relevant question is how to get the maximum increase of permeability for a g iven amount of 
reactive fluid. Numerical and experimental investigations of reactive flows in porous media [pll l33l. l35l [50L |72| suggest 
that there exists an optimum injection rate, which maximizes the permeability gain for a given amount of fluid. If 
the injection rate is relatively small, surface inundation occurs and the increase in permeability is small. On the other 
hand, for very large injection rates, the reactant is exhausted on a uniform opening of the fracture, which is inefficient 
in terms of permeability increase. The optimum flow rate must give rise to spontaneous channeling, since the reactant 
is then used to create a small number of highly permeable channels, which transport the flow most efficiently. To 
quantify the optimization with respect to Pe and Da, we measured the total volume of reactive fluid, Vi„j, that 
must be injected into the fracture in order to increase the overall permeability, K, by a factor of 20. In the constant 
injection rate case, 



Vr, 



QTt = PeWDTr 



(50) 



where T is the time needed for the given permeability increase, measured in units of r. Inserting the definition of r 
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FIG. 15: Permeability as a function of time during dissolution at constant flow rate; results are shown for Pe — 8 and 32. The 
lines correspond to: Da = 1/160 (dotted), Da = 1/40 (dash-and-dotted), Da — 1/10 (dashed), and Da = oo (solid). 

1(15)1 gives 

V mj = TPe^Vo, (51) 

where 

V Q = WLh^^ (52) 

is the volume of reactive fluid needed to dissolve a solid volume equal to the initial pore space in the fracture. Fig. 1161 
shows contour plots of Vi n j/Vo in the Pe — Da plane. A comparison with the dissolution patterns in Fig. 2] suggests 
that optimal injection rates (Pe w 10 — 100, Da > 1) do indeed correspond to a regime of strong channeling. 

The values of Pe and Da cannot be varied independently in the same system, since both the diffusion constant 
and reaction rate are material properties. Changing the injection rate moves the system along a line of constant 
PeDa = kho/D, as shown by the dashed line in Fig. (TBI This produces a characteristic U-shaped dependence of Vt n j 
on Damkohler number, an example of which is shown in the inset to Fig. 1161 The minimum in this curve corresponds 
to the optimal injection rate for a g iven value of PeDa. A similar dependence of Vi n j on Damkohler number has been 
reported previously 0, H (H Hfl Izl • 

An important practical observation is that the optimum flow rate in the mass-transfer-limited regime apparently 
occurs at a constant value of the effective Damkohler number (35|. This result, based on Darcy-scale simulations 
of dissolution in two-dimensional porous media is consistent with the results in Fig. [T|)J When Da > 1, the line of 
constant Da e ff — 1/10 (from Fig.[8j runs near the valley of minimal Vj n j(Pe, Da), However, when the reaction rate 
is reduced, the optimum path shifts to higher Da e tf, as can be seen in the inset to Fig. 1161 



IX. WORMHOLE COMPETITION AND COARSENING OF THE PATTERN 



In Fig. [T71 the dissolution patterns for a larger fracture (800 x 800 pixels) at Pe — 32 and Da — > oo (constant 
pressure drop) are captured at three different instances, corresponding to an aperture increase Ah = 0.15/io, 0.5hty 
and 2ho respectively. Only a small fraction of the channels present at Ah = O.lbho persist to later times (Ah = 
0.5/io); the channels that do survive have advanced far ahead of the dissolution front. The process of channel 
competition is self-similar, and the characteristic length between active (growing) wormholes increases with time, while 
the number of active channels decreases; these systems have been shown to exhibit nontrivial scaling relations (84j . 
The competition between the emerging fingers leads to hierarchical structures that are characteristic of many unstable 
growth processes [H, UK, [3^, H3, [HI] from viscous fingering (75[ and dendritic side-branches growth in crystallization [l6j] 
to crack propagation in brittle solids [47|] . However, due to the finite size of the fracture system, the competition ends 
as soon as the fingers break through to the outlet. In fact, even before breakthrough the simple hierarchical growth 
pattern is disrupted, as shown in the rightmost panel of Fig. 1171 The large pressure gradient at the tips of the leading 
channels causes them to split into two or more daughter branches (l7l . (33l . (3| . 
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FIG. 16: Contour plot of the volume of reactant Vmj needed for a 20-fold increase in permeability. The contours are normalized 
by Vb, the volume of reactant needed to dissolve a solid volume equal to the initial pore space in the fracture. The dashed 
line corresponds to varying the injection rate at PeDa = 0.2, and the inset shows the cross-section of the Pe — Da' 1 surface 
along this line. The dot-dashed line corresponds to Da e ff = 1/10 and indicates a near optimum injection condition in the 
mass-transfer-limited regime. 




FIG. 17: (Color online) Erosion of the lower surface (initially flat) in an artificial fracture at Pe = 32 and Da — > oo, captured 
at (from left to right) Ah = 0.15/io, 0.5/io and 2ho respectively. 
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FIG. 18: (Color online) Contours of the integrated flow field, v 2d = V v 2d ■ v 2d , (left) and its lateral component, Vy d , (right) 
in a small section of the fracture. In the left panel, the lightest/red shading indicates the regions of highest flow, followed by 
green, blue, and black. In the right panel the shading corresponds to the direction of the lateral fluid velocity, white for Vy d > 
and black for Vy d < 0. The integrated flow field, v M , is defined in Eq. (f22|) . 




x 1 



FIG. 19: Two dissolution channels in the fracture (left panel) and the corresponding pressure drops (right panel). The flow 
lines are converging towards a larger channel at the inlet and diverging near the tip of the conduit. (This figure was taken from 
Szymczak and Ladd HU). 

The interaction between wormholes, which underlies the selection process, can be investigated by analyzing the 
flow patterns in the dissolving fracture. Figure [18] shows a magnified view of a part of the sample containing just 
a few channels. The left panel shows the magnitude of the fluid flow in the system, v 2d — V v 2d ■ v 2d ; the flow is 
focused in a few active channels while the rest of the medium is bypassed. The right panel shows the lateral (Vy d ) 
component of the flow: in the white regions v 2d > whereas in the black regions Vy d < 0. It is observed that the 
longest channel, positioned in the center of the sample, is draining flow from the two channels above it as well as from 
the one underneath; v 2d > below the channel and v 2d < above it. Similarly, the second largest channel, situated 
in the lower part of the fracture is draining flow from the smaller channel immediately above it. Thus, out of the 
four large channels in Fig. [HI only two are really active; the other two just supply flow to the active channels. A 
careful examination of the right panel of Fig. [TSl shows that at the downstream end of the channels the flow pattern is 
reversed; the flow is now diverging from the active channels rather than converging towards them as in the upstream 
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FIG. 20: Left: A cross-section of the flow map v1 d (y) at Ah = ho (solid) and Ah = 2h (dashed) (Pe = 8 and Da = 0.1). 
Only four channels have increased their flow rates significantly in the corresponding time interval and at Ah = 2ho about 55% 
of the total flow is focused in the main channel. Center: The fraction of flow focused in the main channel as a function of 
the dissolved volume Ah /ho. The change of slope marks the point when the channel breaks through to the other side of the 
fracture. Right: The ratio of the flow focused in two neighboring channels (A and B) as a function of Ah 



part. 

Fluid flow in the vicinity of the channels is shown schematically in the left panel of Fig. [T9l and can be understood 
from the pressure drop in the channels [84[ , sketched in the right panel. With a constant total pressure drop between 
the inlet and outlet, the pressure gradient in the long channel is steeper than in the short channel, because the flow 
rate is higher. In the upstream part of the fracture the short channel is therefore at a higher pressure than the long 
one, so flow in the surrounding matrix is directed towards the long channel. Downstream the situation is reversed; the 
region around the tip of the long channel is at a higher pressure than the surrounding medium and so flow is directed 
away from the channel, resulting in the converging-diverging flow pattern seen in the simulations. The larger the 
difference in channel lengths, the higher the pressure difference between the channels. The diverging flow pattern at 
the tips leads to the characteristic splitting seen near the breakthrough region in Fig. [T7l where the pressure gradients 
at the tips are large. 

The higher mass flow rates in the longer channels lead to more rapid dissolution and this positive feedback causes 
rapid growth of the longer channels and starvation of the shorter ones. Channel competition is explicitly illustrated in 
Fig- HOI which compares channel lengths and flow rates at different stages of the dissolution process. At the beginning, 
about 20 channels were spontaneously formed by the initial instability at the dissolution front, while at Ah sa h Q (left 
panel) only four of them remain active. As dissolution progresses these remaining four channels self select, and at 
Ah = 2hg a single active channel transports more than a half of the total flow through the sample (center panel). It 
can be seen that channel A which starts out only slightly longer than channel B, drains flow from channel B (right 
panel), slowly at first, but eventually the volumetric flow in channel A becomes an order of magnitude larger than 
that in B. This process repeats itself until only a single channel remains or breakthrough occurs. 



X. INFLUENCE OF THE INITIAL TOPOGRAPHY ON WORMHOLE FORMATION 



It was shown in Sees. IIVI and I Villi that initial topographies characterized by a large number of contact points 
(C = 0.5), relatively high roughness, and small mean aperture, have a well-defined range of Peclet and Damkohler 
numbers where spontaneous channeling is strong. However, in the other geometric extreme, when fracture surfaces 
are far from each other and the roughness is relatively small, the initial inhomogeneities in flow and aperture are 
smoothed out as dissolution proceeds. The effect of fracture roughness can be illustrated by comparing dissolution 
patterns computed in the original fracture geometry, with those obtained after introducing an additional separation 
(extra aperture) h between the surfaces, while keeping the geometry (£ = 0.5), Peclet and Damkohler numbers the 
same. In this case the obstacles do not span the whole aperture, and the relative roughness is a function of both ( 
and h, 



1/2 



/(C, h) 



(l-Qho + h 



(53) 



Figure I2T1 shows the influence of additional aperture on the dissolution patterns for Pe = 8 and Pe = 32, where 
the channeling was found to be strongest. We use an identical arrangement of obstacles in each case (£ = 0.5), and 
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FIG. 21: (Color online) Dissolution patterns for the original topography (a) and for fractures with additional separation between 
the surfaces: h = O.l/io (b) h = 0.5fto (c), h = ho (d), and h = 2ho (e). The flow rates (at constant pressure drop) correspond 
to initial Peclet numbers Pe — 8 (upper panel) and Pe — 32 (lower panel); Da — * oo in both cases. The erosion patterns of an 
initially flat surface were captured at Ah — ho- 



the Damkohler number Da — > oo. The results show a rather well-defined threshold of additional aperture, h tr w ft, , 
beyond which wormholing is strongly suppressed. However, below the threshold the degree of channeling is scarcely 
affected, although the relative roughness is reduced from / = latft = 0to/ = 0.33 at h = ho- It was argued in 
Hanna and Rajaram [4l| that, at larger values of /, flow focusing and channel competition are strongly enhanced, 
whereas for small relative roughness the flow is more diffuse and dissolution becomes uniform. Our results indicate 
a qualitative connection between the roughness of the aperture and channel formation, but not a strong quantitative 
correlation. Despite the large reduction in statistical roughness the dissolution patterns are essentially unchanged up 
to the threshold values of h. In contrast, for separations above the threshold, the patterns change dramatically, with 
the channels disappearing as the value of h is increased from approximately ho to 2/io, yet the statistical roughness 
decreases only slightly, from / = 0.33 to / = 0.2. The absence of a good quantitative correlation between / and the 
degree of channel formation, coupled with the sharp transition to channel suppression suggests that other geometric 
factors, beyond statistical roughness, may play a significant role in determining wormholing. One such geometric 
factor is the degree of contact between the fracture surfaces, which will be examined below. 

The evolution of permeability, Fig. 1221 shows that a small amount of additional separation between the surfaces 
(h = O.lho) does not suppress wormholing (Pe = 32) and can even speed up both dissolution and channel formation 
(Pe = 8). This illustrates two important points. First, a small gap between the fracture surfaces blocks the flow 
pathways almost as effectively as complete contact. Second, a small additional separation enhances the flow rates in 
the spontaneously formed channels, which leads to faster growth, despite the reduction in roughness. 

Further insight into the influence of geometry on the dissolution patterns may be gained by reducing the initial 
coverage, Cj randomly removing some of the protrusions from the original fracture, as illustrated in Fig. 1231 The 
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FIG. 22: Time evolution of the permeability at Pe = 8 (left) and Pe — 32 (right). Results are shown for different additional 
apertures: h — (dotted), h = 0.1/io (dot-dashed), h = 0.5/io (dashed), h — ho (long dashes) and h = 2ho (solid). In the left 
panel the plots for h — ho and h = 2ho overlap and lie almost along the horizontal axis, while in the right panel the plots for 
h = and h = O.l/io overlap. 




FIG. 23: Initial fracture geometries at coverages £ = 0.5 (left), £ = 0.25 (center), and f = 0.05 (right). 



roughness is again reduced, but in a geometrically different fashion from that represented in Fig. [22] Figure l24l 
shows the resulting dissolution patterns at Pe = 32, Da — > oo. Again, there seems to be a well defined threshold 
where channeling is suppressed, occurring at coverages between £ = 0.05 and C = 0.1, which, according to Eq. (|53|) 
corresponds to a relative roughness / w 0.2 — 0.3. Above the threshold, both the diameter of the channels and 
the spacing between them are weakly dependent on the geometry of the system. This supports the notion that 
the characteristics of wormhole formation are primarily functions of Pe and Da, and essentially independent of the 
correlation length of the fracture topography. 



XI. SUMMARY 



In this paper we have studied channeling instabilities in a single fracture, using a fully three-dimensional, microscopic 
numerical method. Channeling was found to be strongest for large reaction rates (mass-transfer-limited regime) and 
intermediate Peclet numbers. We analyzed the observed patterns in terms of two simple models; a Darcy-scale model 
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FIG. 24: (Color online) Dissolution patterns for different coverages (a) £ = 0.5, (b) £ = 0.25, (c) £ = 0.15, (d) £ = 0.05, and 
(e) £ = 0.025. The erosion of an initially flat surface at Pe = 32, Da — > oo, was captured at Aft = 1.25/io. 

for the reaction front and a convection-diffusion model for mass transport in the wormhole. These models qualitatively 
explain the wide variety of dissolution patterns observed in the simulations, and we found quantitative agreement 
in the prediction of the wormhole diameter in the mass-transfer-limited regime. We summarized our simulations in 
terms of a phase diagram separating the different regimes of erosion, and compared our conclusions to experiments 
and other numerical simulations. We also examined the efficiency with which permeability can be increased by acid 
erosion. 

When wormholes form, there is strong competition for the flow, leading to an attrition of the shorter channels. 
We have previously explained a number of detailed observations by a simple network work of the flow in the fracture 
system [84j |. In this paper we have provide new simulation data showing explicitly how the fluid flow is drained 
from the matrix surrounding the dominant channels. This provides confirmation at the microscopic level for key 
assumptions underpinning the network model. 

We found evidence for the existence of a well-defined threshold value of the fracture roughness, / » 0.2 — 0.3, needed 
to trigger a channeling instability. Below that value, the topographic perturbations are smeared out by dissolution 
on similar or shorter timescales than the propagation of the front, and channels do not develop. Above the roughness 
threshold channels do develop, but their size and spacing are controlled by the values of Peclet and Damkohler number, 
and not by the fracture topography. 

Our results are limited to geometries characterized by short-range spatial correlations and lacking the self-affinc 
properties of natural fractures. The analysis of wormholing in such geometries will be the subject of future study. 
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